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Abstract 


Constructing  maps  of  pollution  levels  is  vital  for  air  quality  management,  and  presents 
statistical  problems  typical  of  many  environmental  and  spatial  applications.  Ideally,  such 
maps  would  be  based  on  a  dense  network  of  monitoring  stations,  but  this  does  not  exist. 
Instead,  there  are  two  main  sources  of  information  in  the  U.S.:  one  is  pollution  measurements 
at  a  sparse  set  of  about  50  monitoring  stations  called  CASTNet,  and  the  other  is  pollution 
emissions  data.  The  pollution  emissions  data  do  not  give  direct  information  about  pollution 
levels,  but  instead  are  combined  with  numerical  models  of  weather  and  the  emissions  process 
and  information  about  land  use  and  cover  (collectively  called  Models-3),  to  produce  maps. 

Here  we  develop  a  formal  method  for  combining  these  two  sources  of  information.  We 
specify  a  simple  model  for  both  the  Models-3  output  and  the  CASTNet  observations  in  terms 
of  the  unobserved  ground  truth,  and  estimate  the  model  in  a  Bayesian  way.  This  yields  solu¬ 
tions  to  the  spatial  prediction,  model  validation  and  bias  removal  problems  simultaneously. 
It  provides  improved  spatial  prediction  via  the  posterior  distribution  of  the  ground  truth, 
allows  us  to  validate  Models-3  via  the  posterior  predictive  distribution  of  the  CASTNet  ob¬ 
servations,  and  enables  us  to  remove  the  bias  in  the  Models-3  output  by  estimating  additive 
and  multiplicative  bias  parameters  in  the  model.  We  apply  our  methods  to  data  on  SO2 
concentrations. 

Key  words:  air  pollution,  Bayesian  inference,  change  of  support,  likelihood  approaches, 
Matern  covariance,  nonstationary  process,  spatial-temporal  statistics. 
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1  Introduction 


Emission  reductions  were  mandated  in  the  Clean  Air  Act  Amendments  of  1990  with  the  ex¬ 
pectation  that  they  would  result  in  major  reductions  in  the  concentrations  of  atmospherically 
transported  pollutants.  Maps  of  loadings  of  pollutants  to  aquatic  and  terrestrial  ecosystems 
are  useful  over  different  geo-political  boundaries,  to  discover  when,  where,  and  to  what  extent 
the  pollution  load  is  improving  or  declining.  Ideally,  such  maps  would  be  based  on  a  dense 
network  of  monitoring  stations,  covering  most  of  the  U.S.,  at  which  fluxes  and  concentra¬ 
tions  of  air  pollutants  would  be  measured  on  a  regular  basis.  Unfortunately,  such  a  network 
does  not  exist.  Instead,  there  are  two  main  sources  of  information  about  pollution  levels 
in  the  U.S.,  and  two  resulting  ways  of  constructing  pollution  maps.  The  hrst  is  a  sparse 
set  of  about  50  irregularly  spaced  sites  in  the  Eastern  U.S.,  the  Clean  Air  Status  Trend 
Network  (CASTNet),  at  which  the  EPA  regularly  measures  concentrations  and  fluxes  of 
different  atmospheric  pollutants  (see  Figure  1).  It  would  be  possible  to  use  an  interpolation 
method  to  produce  a  pollution  map.  However,  the  air  pollutants  fluxes  and  concentrations 
are  functions  of  terrain,  atmospheric  turbulence,  vegetation,  the  rate  of  growth  of  the  vege¬ 
tation,  and  other  soil  and  surface  conditions.  Because  these  factors  vary  abruptly  in  space 
and  time  and  because  the  monitoring  stations  are  too  far  from  each  other,  interpolation  of 
the  CASTNet  monitoring  data  is  recognized  to  be  inadequate  for  the  problem  (Clarke  and 
Edgerton,  1997). 

The  second  source  of  information  is  pollution  emissions  data.  The  point  and  area  sources 
emissions  are  available  from  known  sources  of  pollution  such  as  chemical  plants,  generally 
in  the  form  of  annual  totals.  If  the  emissions  data  were  accurate  and  available  at  a  hne 
time  resolution,  and  if  we  had  precise  information  about  local  weather,  land  use  and  cover, 
and  pollution  transport  dynamics,  we  could  in  principle  work  out  pollution  levels  at  each 
point  in  time  and  space  quite  accurately.  This  ideal  is  far  from  being  attained.  However, 
the  available  emissions  data  have  been  combined  with  numerical  models  of  local  weather 
(the  Mesoscale  Model  version  5  (MM5)),  the  emissions  process  (the  Sparse  Matrix  Operator 
Kernel  Emissions  (SMOKE)  model),  as  well  as  information  about  land  use  and  cover,  to 
estimate  pollution  levels  in  space  and  time  (the  Community  Multiscale  Air  Quality  (CMAQ) 
output)  and  to  produce  maps  (Dennis  et  al,  1996).  These  are  not  statistical  models  but  rather 
numerical  deterministic  simulation  models  based  on  systems  of  differential  equations  that 
attempt  to  represent  the  underlying  physics;  they  take  the  form  of  huge  blocks  of  computer 
code.  The  combination  of  these  models  is  referred  to  as  “Models-3”  (models  generation  3). 
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The  models  are  run  by  the  EPA  and  individual  U.S.  states,  and  they  provide  estimates  of 
pollutant  concentrations  and  fluxes  on  regular  grids  in  parts  of  the  U.S.  (see  Figure  2). 

The  output  of  Models-3  generates  averaged  concentrations/fluxes  over  regions  of  dimen¬ 
sions  36km  X  36km.  This  approach  also  may  be  unsatisfactory  for  two  main  reasons.  First, 
the  underlying  emissions  data  are  often  not  of  high  quality  (Dolwick  et  a/.,  2001).  Second, 
the  underlying  models  may  be  inadequate  in  various  ways.  It  seems  clear  that  combining  the 
two  main  approaches  and  sources  of  information,  the  model  estimates  and  the  point  mea¬ 
surements,  could  lead  to  a  better  solution.  So  far,  efforts  to  do  this  have  focused  on  model 
validation,  in  which  model  predictions  are  compared  with  measurements,  and  the  models 
are  revised  and  the  outputs  adjusted  if  discrepancies  are  found  (Dennis  et  a/.,  1990).  The 
hnal  maps  are  still  based  on  the  model  output  alone. 

Model  validation  is  tricky  in  this  case,  because  the  model  predictions  and  the  observations 
do  not  refer  to  the  same  spatial  locations,  and  indeed  are  on  different  spatial  scales.  The 
fact  that  they  are  on  different  spatial  scales  is  called  the  “change  of  support”  problem.  The 
model  predictions  are  averages  over  large  36km  x  36km  grid  squares,  while  the  observations 
are  at  points  in  space;  the  two  are  thus  not  directly  comparable.  One  approach  to  making 
them  comparable  is  to  apply  interpolation  and  extrapolation  methods  to  the  CASTNet  point 
measurements  so  as  to  produce  empirical  estimates  of  grid  square  averages,  and  then  compare 
those  to  the  model  predictions  (Sampson  and  Guttorp,  1998).  One  difhculty  with  this  is  that 
the  interpolated  grid  square  averages  can  be  poor  because  of  the  sparseness  of  the  CASTNet 
network,  and  so  treating  them  as  grand  truth  for  model  validation  may  be  questionable. 

A  related  problem  is  that  the  comparison  does  not  take  into  account  the  uncertainty  in 
the  interpolated  values.  In  this  paper,  we  develop  a  new  approach  to  the  model  validation 
problem,  and  show  how  it  can  also  be  used  to  remove  the  bias  in  model  output,  and  to 
produce  improved  maps  that  combine  model  predictions  with  observations  in  a  coherent 
way.  We  specify  a  simple  model  for  both  Models-3  predictions  and  CASTNet  observations 
in  terms  of  the  unobserved  ground  truth,  and  estimate  it  in  a  Bayesian  way.  Solutions  to  all 
the  problems  considered  here  follow  directly.  Model  validation  then  consists  of  comparing 
the  CASTNet  observations  with  their  predictive  distributions  given  the  Models-3  output. 
Bias  removal  follows  from  estimation  of  the  bias  parameters  in  the  model.  Maps  of  pollution 
levels  and  of  the  uncertainty  about  them  taking  into  account  all  the  available  information  are 
based  directly  on  the  posterior  distribution  of  the  (unobserved)  ground  truth.  The  resulting 
approach  takes  account  of  and  estimates  the  bias  in  the  atmospheric  models,  the  lack  of 
stationarity  in  the  data,  the  ways  in  which  spatial  structure  and  dependence  change  with 
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Models-3:  S02  Concentrations 


Figure  2:  Output  of  Models-3,  weekly  average  of  SO2  concentrations  (ppb),  for  the  week  of 
July  11,  1995.  The  resolution  is  36  k-m?. 
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locations,  the  change  of  support  problem,  and  the  uncertainty  about  these  factors.  It  is 
an  instance  of  the  Bayesian  melding  framework  for  inference  about  deterministic  simulation 
models  (Poole  and  Raftery,  2000),  and  its  implementation  is  quite  straightforward. 

In  Section  2  we  describe  the  statistical  model,  in  Section  3  we  show  how  to  estimate  it, 
and  in  Section  4  we  show  some  of  our  results  for  model  validation  and  map  construction 
using  the  combined  data  for  the  air  pollution  problem. 

2  The  Statistical  Model 

2.1  General  Framework 

We  do  not  treat  CASTNet  measurements  as  the  “ground  truth” .  Instead,  we  assume  that 
there  is  some  smooth  underlying  (but  unobserved)  held  that  measures  the  “true”  concen- 
tration/hux  of  the  pollutant  at  each  location.  CASTNet  data  are  these  “true  values”  plus 
some  measurement  error.  The  Models-3  output  can  also  be  written  in  terms  of  this  true  un¬ 
derlying  (unobservable)  process,  with  some  parameters  that  explain  the  bias  and  microscale 
noise  in  Models-3.  The  truth  is  assumed  to  be  a  smooth  underlying  spatial  process  with 
some  parameters  that  explain  the  large  scale  and  short  scale  dependency  structure  of  the  air 
pollutants. 

Our  objectives  are  model  validation  and  bias  removal  for  Models-3,  and  construction  of 
reliable  maps  of  air  pollution  combining  Models-3  and  CASTNet.  We  validate  Models-3 
by  obtaining  the  posterior  predictive  distribution  of  CASTNet  given  Models-3  output.  We 
remove  the  bias  in  Models-3  by  obtaining  the  posterior  distribution  of  the  bias  parameters 
given  CASTNet  data  and  Models-3  output.  We  construct  reliable  maps  of  air  pollutants 
simulating  values  from  the  posterior  predictive  distribution  of  the  true  values  (underlying 
process)  given  CASTNet  data  and  Models-3  output. 

2.2  Statistical  Models  for  CASTNet  and  Models-3  Output 

Our  general  modeling  framework  is  shown  in  Figure  3.  We  do  not  consider  CASTNet 
measurements  to  be  the  “ground  truth”,  because  there  is  measurement  error.  Thus,  we 
assume  there  is  an  underlying  (unobserved)  held  Z(s),  where  Z(s)  measures  the  “true” 
concentration/hux  of  the  pollutant  at  location  s.  At  station  s  we  make  an  observation  .^(s), 
corresponding  to  the  CASTNet  observation  at  this  station,  and  we  assume  that 

Z(s)  =  Z(s) +  e(s),  (1) 
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Figure  3:  General  Modeling  Framework. 

where  e(s)  ~  iV(0,ag)  represents  the  measurement  error  (nugget)  at  location  s.  The  process 
e(s)  is  independent  of  Z(s). 

The  true  underlying  process  Z  is  a  spatial  process  with  a  nonstationary  covariance, 

Z{s)  =/r{s)  +  e(s),  (2) 

where  Z(s)  has  a  spatial  trend,  /r(s),  that  is  a  polynomial  function  of  s  with  coefhcients  /3.  We 
assume  that  Z{s)  has  zero-mean  correlated  errors  e(s).  The  process  e(s)  has  a  nonstationary 
covariance  with  parameter  vector  6  that  might  change  with  location. 

We  could  model  the  output  of  the  EPA  physical  models  as  follows: 

Z{s)  =  a{s)  +  b{s)Z{s)  +  6{s).  (3) 

Here,  the  parameter  function  a(s)  measures  the  additive  bias  of  the  air  quality  models  at 
location  s,  and  the  parameter  function  b{s)  accounts  for  the  multiplicative  bias  in  the  air 
quality  models.  The  process  5{s)  ~  A^(0,a|)  explains  the  random  deviation  at  location  s 


6 


with  respect  to  the  underlying  true  process  Z(s).  The  process  S(s)  is  independent  of  Z(s) 
and  e(s),  which  is  the  error  term  for  CASTNet.  Since  the  outputs  of  Models-3  are  not  point 
measurements  but  areal  estimations  in  subregions  Bi, that  cover  the  domain,  D,  we 
have 

Z{Bi)  =  f  a{s)ds  +  b  f  Z{s)ds+  f  5{s)ds  (4) 

J  Bi  J  Bi  J  Bi 

for  i  =  1, . . . ,  m.  We  model  the  function  a(s)  as  a  polynomial  in  s  with  a  vector  of  coefhcients, 
oo,  and  6  is  a  unknown  constant  term. 

For  spatial  prediction  we  simulate  values  of  Z  from  the  posterior  predictive  distribution: 

P{Z\Z,Z).  (5) 

For  model  validation  we  simulate  values  of  CASTNet  given  models-3,  from  the  following 
posterior  predictive  distribution: 

P{Z\Z,a  =  0,b  =  1).  (6) 

For  bias  removal  we  simulate  values  of  the  parameters  a  and  b  from  the  posterior  distri¬ 
bution: 


P{a,b\Z,Z). 


(7) 


2.3  Change  of  Support 

The  change-of-support  problem  occurs  when  we  combine  data  sources  with  different  sup¬ 
ports,  or  when  the  supports  of  predictand  and  data  are  not  the  same.  Here,  we  have  point 
measurements  at  the  CASTNet  sites,  and  then  we  observe  the  output  of  Models-3  averaged 
over  regions,  Hi, ... ,  of  dimensions  36km  x  36km. 

In  this  section  we  discuss  algorithms  to  calculate  the  covariance  for  areal  measurements 
and  the  posterior  predictive  distribution  of  a  random  process  at  a  point  location  .Z'(xo)  given 
data  on  block  averages,  Z{Bi), . . . ,  Z{Bm),  where  some  of  the  blocks  might  be  just  a  point. 
The  covariance  for  the  block  averages  is 

C0Y{Z{Bi),Z{Bj))  =  f  f  C{u,v)dudv/\Bi\\Bj\,  (8) 

J  Bi  J  Bj 

where 

C(u,v)  =  cov(Z(u),Z(v)), 
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C  being  a  possibly  nonstationary  spatial  function.  If  Bi  =  Sj  (a  point)  the  covariance  is 
defined  by 

co\{Z{si),Z{Bj))  =  f  C{si,v)dy/\Bj\.  (9) 

J  Bj 

In  practice,  for  each  pixel  Bj  we  draw  an  independent  set  of  locations  Uj^,  r  =  1,  2, . . . , 
uniformly  over  Bj,  and  we  approximate  the  integral  in  (9)  with  the  following  expression: 

cm{Z{B,),Z{B,))=L;^  C{Uir,Ujr').  (10) 

r  r' 

We  use  (9),  approximated  by  (10),  to  derive  the  covariances  of  the  block  averages  Z{Bi), 
i  =  1, . . .  ,iV,  in  terms  of  the  pointwise  covariance  C{u,v).  This  is  then  used  to  dehne  a 
likelihood  function  for  the  parameters  of  the  covariance  function  for  the  process  Z  in  terms 
of  the  observed  block  averages  Z{Bi),  Z{B2), . . . ,  Z{Bj^). 

2.4  Methods  for  Combining  Data  with  Different  Spatial  Resolu¬ 
tions 

The  spatial  processes  Z  and  Z  are  measuring/estimating  the  same  quantity  of  interest  but 
they  have  different  spatial  resolution.  The  process  Z  is  observed  at  n  locations  Si, . . .  ,s„, 
so  we  have  n  point  measurements  Z  =  ^Z(si), . . . ,  Z{sn)^  .  The  process  Z  is  averaged  over 

the  subgrids  Bi, . . . ,  B^,  so  that  we  observe  Z  =  (^Z{Bi), . . . ,  Z{Bfnfj  ■  We  model  Z  and 
Z  in  terms  of  an  underlying  unobservable  process  Z,  which  is  the  true  value  of  the  quantity 
of  interest.  We  take  into  account  the  measurement  error  and  potential  bias  of  Z  and  Z  as 
approximations  to  Z.  Here  we  consider  model  (1)  for  Z  and  model  (4)  for  the  process  Z. 

We  now  deduce  the  joint  distribution  of  Z  and  Z  conditioning  on  the  value  of  the  param¬ 
eters  in  models  (1)  and  (4).  We  could  write  this  distribution  as  a  function  of  the  parameters 
to  calculate  the  MLE  for  the  parameters  in  models  (1)  and  (4).  Since  in  practice  this  cal¬ 
culation  will  be  hard,  we  also  present  a  Bayesian  approach  to  estimate  the  parameters.  We 
have 


where 


~  A/" 


A 

d  +  bfi 


S 

S 


A=  (/^(Sl),...,/^(Sn))^, 


(11) 
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and 


a 


/  a{s)ds, .  . 

. ,  /  a{s)ds  j 

J  Bi 

^  Bm  J 

/  H{s)ds, . . 

JBi 

. ,  /  H{s)ds] 

JBm  / 

In  (11)  S  is  the  covariance  of  Z  which  is  the  covariance  of  Z  plus  a  nugget  effect,  S  is  the 
covariance  of  Z,  and  S  is  the  cross-covariance  between  the  point  measurements  Z  and  the 
block  averages  Z.  We  write  S  to  denote  the  covariance  of  the  multivariate  normal  distribution 
of  ^Z,  Z^  ,  so  that  S  is  a  (n  -I-  m)  X  (n  -I-  m)  matrix  given  by 


'^n+j,i  —  '^i,n+j  —  COV 


V  ) ,  ^(822)^  ^*(^21:^22)  l{ii=*2}W 

=b‘^  C{s,,v,0)dv/|5,|, 

/«  „<^{u,v,6>)dudv 


-‘n+ji  ,n+j2 


=  cov  (^Z{BjJ,Z{Bj^)^  = 


^jl  ^32 


for  ii,  12  <  n, 

for  i  =  1, . . . ,  n  and  j  =  1, . . . ,  m. 

for  jij2  =  1,  •  •  ■ , 


where  C{u^  v,  0)  is  the  covariance  between  Z{u)  and  Z{v)^  and  0  is  the  parameter  associated 
with  the  covariance  of  the  true  underlying  process  Z,  the  function  1a  (x)  is  an  indicator  for 
the  set  yl,  i.e.  takes  the  value  1  when  x  e  A  and  0  otherwise. 

The  goal  is  to  predict  the  value  of  Z  at  location  xq  given  the  data.  Thus  we  need  the 
conditional  distribution  of  .Z'(xo)  given  the  observations,  assuming  all  the  parameters  are 
known.  We  use  the  following  classical  result  of  multivariate  analysis  (see,  e.g.,  Mardia,  Kent 
and  Bibby  (1979),  p.  63).  If  we  consider  a  partitioned  multivariate  normal  vector 


A2 


A7 


{{ 


/^i 

/^2 


Sii  S12 

^21  S22 


)}■ 


(12) 


then  the  conditional  distribution  of  Xi  given  X2  is  normal  with  mean 

Hi  -|-  'Ei2'Z22  iX2  —  H2) 

and  variance 


^12^22^^21. 


Applying  this  result  with  Xi  =  .Z'(xo),  and  A2  =  Z  =  (Z,  Z)^,  then  in  our  previous  notation 
Ell  =  cov(2'(xo),  Z{-xo))  =  CTq,  E22  =  E,  and  E12  =  cov(2'(xo),  Z)  =  r,  where  r  is  a  (n  -|-  m) 
dimensional  vector  with  components, 

n  =  cov{Z(xo),  Zi}  =  cov  |.Z'(xo),  i'(sj)|  =  C'(xo,Sj,0),  for  i  =  1, . . . ,  n. 


m, 
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for  j  =  1, . . . ,  m, 


Tn+j  =  cov{Z{xo),  In+j]  =  cov  |z(xo),  1(5^) |  =  j  C ,d)d-v /\Bjl 

and  Zj  denotes  the  component  of  Z.  We  deduce  then  that  the  conditional  distribution  of 
Z(xo)  given  {Z,  Z}  is  normal  with  mean 

/r(xo)  +  r^S“^(Z  -  /r)^, 


where 


and  variance 


11=  (/t,  a  +  hjl) 


T 


T’sn-l 


ao  -  r  S 


r. 


When  the  goal  is  to  predict  Z  at  a  location  xq,  the  Bayesian  solution  is  the  predictive 
distribution  of  .Z'(xo)  given  the  observations  Z, 

p(Z(xo)|Z)oc  yp{Z(xo)iZ,  0)p(0|Z)(i0.  (13) 

A  Gibbs  sampling  approach  is  used  to  simulate  m  values  from  the  posterior  of  the  vector  pa¬ 
rameter  0,  where  0  =  (ag,  oq,  6,  0).  Thus,  the  predictive  distribution  is  approximated 

by  the  Rao-Blackwellized  estimator. 


p{Z{xo)iZ)  =  —  ^p(Z(xo)iZ,0W), 

TTh 


(14) 


where  is  the  *-th  draw  from  the  posterior  distribution.  We  propose  a  Gibbs  sampling 
approach,  with  three  stages.  In  stage  1  we  sample  from  the  conditional  posterior  of  the 
parameters  {/3,  0),  the  parameters  that  represent  the  lack  of  stationary  of  Z.  In  the  second 
stage  we  sample  from  the  conditional  posterior  of  the  parameters  that  explain  the  bias  of  Z 
and  the  measurement  error  of  Z  and  Z,  namely  (ag,ao,6,  a|).  In  stage  three,  we  simulate 
values  of  Z  (the  unobserved  true  values)  at  the  locations  where  we  have  measurements  from 
the  process  Z  and  Z.  We  cycle  through  the  three  stages. 


2.5  Modeling  a  Nonstationary  Covariance 

The  spatial  patterns  shown  by  the  air  pollutant  fluxes  and  concentrations  change  with  loca¬ 
tion,  so  that  the  underlying  process  Z  in  (2)  is  nonstationary  and  the  standard  methods  of 
spatial  modeling  and  interpolation  are  inadequate.  In  this  section  we  give  a  new  method¬ 
ology  for  spatial  modeling  of  nonstationary  covariance.  More  specihcally,  we  represent  the 
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process  locally  as  a  stationary  isotropic  random  field  with  some  parameters  that  describe 
the  local  spatial  structure  (Fuentes,  2001a, b).  These  parameters  are  allowed  to  vary  across 
space  and  reflect  the  lack  of  stationarity  of  the  process. 

Consider  a  Gaussian  spatial  process  Z(x),  where  x  varies  over  a  domain  D  contained  in 
a  d-dimensional  Euclidean  space  for  some  d  >  1.  Typically,  d  =  2.  We  represent  Z  as  a 
convolution  of  local  stationary  processes  (Fuentes  and  Smith,  2001): 

Z(x)  =  f  iF(x  -  s)Z0(s)(x)ds,  (15) 

Jd 

where  K  is  a  kernel  function  and  x  e  H  is  a  family  of  (independent)  stationary 

Gaussian  processes  indexed  by  6.  The  parameter  d  is  allowed  to  vary  across  space  to  reflect 
the  lack  of  stationary  of  the  process.  The  stochastic  integral  (15)  is  dehned  as  a  limit  (in 
mean  square)  of  approximating  sums  (e.g.,  Cressie,  1993,  p.  107,  Yaglom,  1962,  p.  23). 
Each  stationary  process  Zg(^s)(x)  has  a  mean  function  /Xg  that  is  constant,  i.e.  jUs  does  not 
depend  on  x.  We  propose  a  parametric  model  for  the  mean  of  Z, 

E{Z{x)}  =  ii{x;/3), 

where  ji  could  be  a  polynomial  function  of  x  with  coefhcients  /9. 

The  covariance  of  Ze(s)  is  stationary  with  parameter  0(s), 

COv{Z0(s)(si),  .^0(s)(s2)}  =  C'0(s)(si  —  S2). 


The  process  Ze{s)  could  have  a  Matern  stationary  covariance  (Matern,  1960): 

X  2:.  (2^^y^|x|M)^^/C^,(2t^y^|x|/p,),  (16) 

where  is  a  modihed  Bessel  function  and  d{s)  =  (z/^,  cr^,  p^).  The  parameter  pg  measures 
how  the  correlation  decays  with  distance;  generally  this  parameter  is  called  the  range.  The 
parameter  Ug  is  the  variance  of  the  random  held,  i.e.  Ug  =  var(Z0(s)(x)),  where  the  covariance 
parameter  Ug  is  usually  refereed  to  as  the  sill.  The  parameter  Vg  measures  the  degree  of 
smoothness  of  the  process  Ze(^g,y  The  higher  the  value  of  Vg  the  smoother  would  be; 
e.g.  when  Vg  =  |,  we  get  the  exponential  covariance  function.  In  the  limit  as  Vg  ^  00  we 
get  the  Gaussian  covariance 

G0(s)(x)  = 

The  covariance  G(si,  S2;  9)  of  Z  is  a  convolution  of  the  local  covariances  Ce(s){si  —  S2), 

G(si,S2;0)=  [  K{si  -  s)K{s2  -  s)Ce{s){si  -  S2)ds.  (17) 

Jd 
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In  (17)  every  entry  requires  an  integration.  Since  each  such  integration  is  actually  an 
expectation  with  respect  to  a  uniform  distribution,  we  propose  Monte  Carlo  integration.  We 
propose  to  draw  a  systematic  sample  of  locations  s^,  m  =  1,2,  ...,M  over  D.  Hence,  we 
replace  (7(81,82;  0)  with 

M 

C'm(Si,  S2;  d)  =  M~^  ^  K(8i  -  Sm)K{s2  -  S„)(70(s^) (8i  -  82).  (18) 

m=l 

This  is  a  Monte  Carlo  integration  which  can  be  made  arbitrarily  accurate  and  has  nothing  to 
do  with  the  data  Z.  The  sampling  points  s^,  m  =  1,  2, ...,  M,  determine  subregions  of  local 
stationarity  for  the  process  Z.  We  increase  the  value  of  M  until  convergence  is  achieved. 

3  Estimation 

In  this  section  we  explain  how  to  efhciently  implement  our  algorithm  for  spatial  prediction 
combining  observations  with  the  output  from  numerical  models. 

3.1  Algorithm 

1.  Posterior  predictive  values  for  Z 

For  spatial  prediction  the  quantity  of  interest  is  the  predictive  distribution  for  .Z'(xo) 
given  the  observed  values  Z.  We  use  expression  (14)  to  approximate  the  predictive  dis¬ 
tribution  with  a  Rao-Blackwellized  estimator,  conditioning  on  the  posterior  simulated 
values  for  all  the  parameters,  using  for  this  simulation  the  following  Gibbs  algorithm. 

2.  Algorithm  for  Gibbs  sampling 

We  discuss  now  how  to  sample  from  the  posterior  distribution  of  the  parameters.  In  our 
Gibbs  sampling  approach  there  are  three  stages.  We  alternate  between  the  parameters 
that  measure  the  lack  of  stationarity,  (^,  6)  (Stage  1),  the  parameters  that  measure  the 
bias  of  Models-3  and  the  measurement  error  of  CASTNet  (Stage  2),  and  the  unobserved 
true  values  of  Z  at  all  the  CASTNet  sites  and  at  the  blocks  where  we  have  the  Models-3 
output  (Stage  3). 

Gibbs  sampling:  Stage  1. 

We  obtain  the  conditional  posterior  for  the  parameters  that  measure  the  lack  of  sta¬ 
tionarity,  (/3,  0(s)),  conditioning  on  the  values  of  Z  that  are  updated  in  Stage  3.  The 
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posterior  of  {P,  0(s))  will  be  completely  specified  once  we  define  the  priors  for  (/3,  0(s)), 
because  we  have  that 

[Z\P^  6]  is  Gaussian, 

where  the  brackets  [  ]  are  used  here  to  denote  densities. 

Gibbs  sampling:  Stage  2. 

We  obtain  the  conditional  posterior  for  the  parameters  ao,6,  a|  and  that  measure 
the  bias  and  uncertainty  of  Models-3,  and  the  measurement  error  of  CASTNet. 

The  posterior  of  given  the  n  values  of  Z  and  Z  at  the  CASTNet  sites  (updated  in 
Stage  3),  can  be  easily  obtained,  because  we  have  the  following  regression  problem: 

Z(s)  =  Z{s)  +  e(s), 

where  is  the  variance  of  the  error  term  e(s),  and  Z{s)  is  independent  of  e(s).  We 
have  that 


[Z{s)\Z{s)^a^]  is  normal  with  mean  Z{s)  and  variance 
Then,  the  posterior  of  is  proportional  to 

[Z{xi), . . . ,  Z(x„)|Z(xi), . . . ,  Z(x„),  ag][ag] 

where  [al]  denotes  the  prior  distribution  for  ag,  and  xi, . . .  ,x„  are  the  n  CASTNet 
sites. 

The  posterior  distributions  of  ao,b,ag  given  the  values  of  Z  and  Z  (updated  in  Stage 
3)  at  the  m  blocks,  can  be  easily  calculated,  because  we  have  the  following  regression 
problem: 

Z{Bi)  =  /  a{s)ds  +  b  /  Z{s)ds+  /  5{s)ds^ 

J  Bi  J  Bi  J  Bi 

where  is  the  variance  of  the  error  term  5(s),  and  Z{s)  is  independent  of  5(s).  It 
follows  that 

[Z{Bi), . . .  ,Z{Bm)\Z{Bi), . . .  ,Z{Bm),ao,b,aj] 

is  normal  with  mean  a+b  {Z{Bi), . . . ,  Z{Bfn)},  where  a  =  |  a(x)(ix,  ■  ■  ■ ,  a(x)(ix|  , 
and  a  diagonal  covariance  matrix  with  diagonal  elements  o‘l\Bi\.  Thus,  the  posterior 
of  oo,  b,  ag  is  proportional  to 

[Z{Bi), . . .  ,Z{Bm)\Z{Bi), . . .  ,Z{Bm),ao,b,aj][ao,b,aj]. 
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Gibbs  sampling:  Stage  3. 

We  simulate  values  of  Z  (the  unobserved  true  values)  at  the  n  locations  where  we  have 
measurements  for  Z,  and  also  at  the  m  blocks  where  we  observe  Z,  conditioning  on  the 
values  of  ^,0  (updated  in  Stage  1)  and  Z.  The  simulated  values  at  the  m  blocks  are 
obtained  by  simulating  values  of  Z  at  a  sample  of  locations  within  each  block.  Then 
Z{Bj)  is  approximated  by  L~^  where  Sj^, . . . ,  is  a  centered  systematic 

sample  in  Bi. 

For  model  validation,  we  simulate  values  from  the  posterior  distribution  of  CASTNet 
given  Models-3, 

P{Z\Z,  a  =  0,b  =  1), 

and  we  compare  the  actual  observations  with  this  simulated  posterior  predictive  distribution. 

4  Application:  Air  Pollution  Data 

The  regional  scale  air  quality  models  (Models-3)  run  by  the  U.S.  EPA  estimate  hourly 
concentrations  and  fluxes  of  different  air  pollutants.  The  primary  objective  of  Models-3  is  to 
improve  the  environmental  management  community’s  ability  to  evaluate  the  impact  of  air 
quality  management  practices  for  multiple  pollutants  at  multiple  scales,  as  part  of  the  process 
of  regulating  air  pollution.  The  spatial  domain,  D,  is  a  regular  grid  (81x87),  the  dimensions 
of  each  pixel  in  the  grid  are  36km  x  36km.  Models-3  provides  hourly  concentrations  for 
each  pixel.  As  an  example  we  examine  sulfur  dioxide.  Figure  2  shows  the  weekly  averaged 
concentrations  of  SO2  for  the  week  starting  July  11,  1995.  Our  hrst  objective  is  the  validation 
of  Models-3. 

The  Clean  Air  Status  and  Trends  Network  (CASTNet)  measures  SO2  weekly  averaged 
concentrations  and  fluxes  at  50  sites  (see  Figure  1).  We  propose  to  combine  both  sources  of 
information,  namely  the  monitoring  data  and  the  output  of  the  air  quality  models,  to  get 
more  reliable  maps  of  air  pollutants  concentrations  and  fluxes. 

In  Figure  4(a)  we  show  the  CASTNet  values  at  6  selected  sites  that  are  representative 
of  different  metereological,  land  use,  and  altitude  conditions.  For  validation  of  Models- 
3,  in  Figure  4(b)  we  show  the  modes  of  the  posterior  predictive  distribution  (eq.  (6))  of 
CASTNet  given  Models-3  at  the  6  selected  sites.  There  is  clear  evidence  that  Models-3 
is  overestimating  the  concentrations  of  S02-  We  modeled  the  covariance  for  Models-3  using 
equation  (17),  taking  into  account  the  lack  of  stationarity  and  the  change-of-support  problem 
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S02  concentrations  (CASTNet) 


S02  concentrations  (Models-3) 


(b) 


Figure  4:  (a):  Weekly  average  of  SO2  concentrations  (ppb)  at  some  selected  CASTNet  sites, 
for  the  week  of  July  11,  1995.  (b):  Weekly  average  of  SO2  concentrations  (ppb)  from  Models-3 
interpolated  at  6  selected  locations,  for  the  week  of  July  11,  1995,  for  Models-3  validation. 


(we  calculated  the  covariances  involving  block  averages  by  drawing  a  set  of  4  locations  in 
each  pixel). 

Now,  we  use  the  methodology  presented  in  this  paper  to  estimate  the  bias  in  Models-3 
for  bias  removal,  using  expression  (7).  We  modeled  Models-3  in  terms  of  an  underlying 
unobservable  process  Z  with  the  true  values  of  SO21  but  we  added  an  additive  constant  bias, 
a  multiplicative  constant  bias,  and  a  measurement  error  term.  We  also  modeled  CASTNet 
in  terms  of  the  “true”  process  Z  and  we  added  a  measurement  error  term  (see  Section  2.2). 

Figures  5  and  6  show  the  posterior  distributions  of  some  covariance  parameters  for  the 
underlying  process  Z  at  the  selected  sites  shown  in  Figure  4  (a).  We  used  vague  gamma 
priors  for  all  the  Matern  covariance  parameters,  except  for  the  sill  parameter  for  which  we 
used 

p{a)  oc 

which  is  a  uniform  prior  for  log(a).  The  sill  parameter  changes  with  location  as  illustrated 
by  the  variation  in  the  distributions  in  Figure  6.  Thus,  this  indicates  a  lack  of  stationarity. 
The  range  parameter  does  not  change  much  with  location  (Figure  5).  The  smoothing  pa¬ 
rameter  does  not  change  with  location  either,  and  is  always  close  to  1/2  (exponential).  We 
implemented  the  nonstationary  model  (15)  with  weight  function  K{u  —  s)  =  Mo  (^)  , 
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Range  Parameter 


0  50  100  150  200  250 


Posterior  distributions 

Figure  5:  Posterior  distributions  for  the  range  parameter  (km)  of  the  Matern  covariance  for 
the  SO2  concentrations  of  Z,  for  the  week  starting  July  11,  1995,  at  the  6  selected  locations 
showed  in  Figure  4  (a). 

where  -K’o(u)  is  the  quadratic  weight  function 

^o(u)  =  ^(1 -Mi^)+^(1  -  V)+,  (19) 

for  u  =  (mi,  M2).  The  bandwidth  parameter  h  is  dehned  as  //2  +  //2e,  where  I  is  the  distance 
between  the  sample  points  si, . . . ,  sm  in  (18),  and  e  is  a  value  between  0  and  1.  For  e  we  used 
a  uniform  prior  in  the  interval  [0, 1].  The  parameter  e  determines  the  amount  of  overlapping 
between  the  subregions  of  stationarity  centered  at  the  sampling  points  si, . . . ,  sm,  and  h  can 
be  interpreted  as  the  diameter  of  the  subregions  of  stationarity. 

The  mode  of  the  posterior  distribution  for  the  parameter  that  measures  the  measurement 
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Figure  6:  Posterior  distributions  for  the  sill  parameter  of  the  Matern  covariance  for  the  SO2 
concentrations  of  Z,  for  the  week  starting  July  11,  1995,  at  the  6  selected  locations  showed 
in  Figure  4  (a). 
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error  for  CASTNet  is  .8,  and  for  Models-3  it  is  .1.  The  mode  of  the  posterior  distribution 
for  the  parameter  that  measures  the  multiplicative  bias  for  Models-3  is  .5  with  a  standard 
error  of  .5,  and  for  the  additive  bias  we  have  a  polynomial  of  degree  4. 

We  use  a  Bayesian  approach  for  spatial  prediction.  We  sample  from  the  predictive  dis¬ 
tribution  (eq.  (5))  of  the  true  underlying  process  Z  given  the  observations  and  the  model 
output,  taking  into  account  the  lack  of  stationarity  and  the  change-of-support  problem  (we 
calculated  the  block  averages  covariances  drawing  a  set  of  4  locations  in  each  pixel).  We 
use  posterior  predictive  checks  (PPG)  as  suggested  by  Rubin  (1984)  and  Raftery  (1988)  for 
validation  of  the  proposed  model  for  Z.  Thus,  we  compare  the  posterior  predictive  distri¬ 
butions  of  the  “true”  process  Z  at  different  locations  to  the  observed  data,  and  we  judge 
whether  the  generated  data  are  similar  to  the  CASTNet  data  (Figure  4(a)).  Figure  7  shows 
simulated  values  from  the  posterior  predictive  distribution  (eq.  (5))  of  SO2  weekly  average 
concentrations  at  the  CASTNet  selected  sites.  Figure  8  is  a  comparison  of  the  generated 
data  through  the  Bayesian  melding  approach  proposed  in  this  paper  to  CASTNet,  taking 
into  account  the  uncertainties  about  CASTNet  and  Models-3.  The  graph  on  the  right  in 
Figure  8  shows  the  CASTNet  measurements  for  the  week  starting  July  11,  1995,  versus  the 
modes  and  90%  credible  intervals  of  the  predictive  Bayesian  distributions  (eq.  (5))  derived 
from  the  Bayesian  melding  approach  at  the  CASTNet  locations.  The  dotted  lines  indicate 
a  90%  credible  region  for  the  CASTNet  values.  All  the  modes  fall  within  the  credible  region 
for  CASTNet,  so  that  the  generated  data  are  similar  to  the  observed  data.  On  the  other 
hand,  the  graph  on  the  left  in  Figure  8  shows  CASTNet  measurements  versus  the  modes 
and  90%  credible  intervals  of  the  predictive  Bayesian  distributions  (eq.  (6))  derived  only 
from  Models-3  (without  combining  Models-3  with  CASTNet)  at  the  CASTNet  locations  for 
validation  of  Models-3.  The  modes  in  the  latter  plot  do  not  fall  within  the  credible  bands  for 
CASTNet.  This  hgure  shows  the  improvement  obtained  in  the  prediction  of  SO2  by  com¬ 
bining  CASTNet  and  Models-3  through  the  Bayesian  melding  approach  presented  in  this 
paper. 

In  Table  1  we  have  the  modes,  sample  standard  deviations,  and  90%  credible  intervals 
of  the  posterior  predictive  distribution  (5)  for  SO2  at  the  6  selected  sites.  As  expected,  we 
get  very  high  variability  at  the  Indiana  site.  This  site  is  very  close  to  several  coal  power 
plants,  and  so  the  SO2  levels  can  be  very  high  or  very  low  depending  on  wind  speed,  wind 
direction,  and  on  the  atmospheric  stability.  The  sites  in  Maine  and  Florida  have  the  lowest 
SO2  levels  and  variability.  The  agricultural  site  in  Illinois  and  the  site  in  North  Carolina 
have  similar  behavior  in  terms  of  SO2  levels.  The  site  in  North  Carolina  is  not  far  from  the 
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Figure  7:  Predictive  posterior  distributions  for  the  SO2  concentrations  of  Z,  at  the  6  selected 
locations  showed  in  Figure  4,  for  the  week  starting  July  11,  1995.  The  predictive  posterior 
values  are  for  the  spatial  underlying  process  Z  given  CASTNet  and  Models-3. 
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Table  1:  Columns  2-5  in  this  table  show  the  modes,  standard  errors  and  90%  credible  intervals 
of  the  posterior  predictive  distribution  (eq.  (5))  for  the  underlying  process  Z  measuring  the 
true  SO2  concentrations  at  the  6  selected  sites.  Column  6  are  the  CASTNet  values  (Z). 
Columns  7-9  show  the  modes  and  the  corresponding  90%  credible  intervals  of  the  posterior 
predictive  distribution  (eq.  (6))  for  model  validation. 


Site 

Mode 

S.E. 

90% 

C.  I. 

CASTNet 

Models-3 

90% 

C.  I. 

Maine 

0.18 

0.08 

0.15 

0.25 

0.15 

0.33 

0.10 

0.43 

Illinois 

2.80 

0.84 

2.55 

3.41 

3.29 

3.33 

2.17 

5.03 

North  Carolina 

1.38 

0.98 

1.18 

2.08 

0.90 

5.32 

3.67 

6.67 

Indiana 

0.98 

5.31 

0.74 

5.63 

3.14 

9.59 

4.20 

20.50 

Florida 

0.91 

0.16 

0.87 

1.05 

0.57 

0.52 

0.20 

0.80 

Michigan 

0.83 

0.40 

0.79 

1.14 

1.02 

1.04 

0.53 

1.70 

Tennessee  power  plants,  and  the  site  in  Illinois  is  also  relatively  close  to  some  Midwestern 
power  plants.  The  site  in  Michigan,  which  is  very  close  to  Lake  Michigan  and  relatively  far 
from  power  plants,  also  has  low  SO2  levels. 

In  Table  1  we  also  show  the  CASTNet  values,  to  judge  if  the  generated  data  are  similar 
to  the  CASTNet  data.  Considering  that  the  uncertainty  about  CASTNet  is  0.8  ppm,  the 
predictive  values  are  fairly  similar  to  CASTNet,  for  the  most  part.  The  CASTNet  values 
in  table  1  at  the  6  sites  represent  the  5th,  91st,  1st,  68th,  0th,  and  79th  percentiles  of 
the  corresponding  posterior  predictive  distributions  from  the  Bayesian  melding  approach. 
The  last  3  columns  in  Table  1  are  the  modes  and  the  corresponding  90%  credible  intervals 
of  the  posterior  predictive  distribution  (6)  for  model  validation,  which  are  also  plotted  in 
Figure  4  (b).  We  can  clearly  appreciate  the  bias  in  the  numerical  models  by  comparing  these 
values  with  CASTNet.  For  instance,  for  the  site  in  North  Carolina,  the  interpolated  value 
of  Models-3  is  5.32  ppm  (s.e.  3.00  ppm)  while  the  CASTNet  value  is  only  0.90  ppm  (see 
Figure  8).  We  could  remove  the  bias  in  these  interpolated  Models-3  values  by  taking  into 
account  the  additive  bias  measured  by  a(x)  (a  polynomial  of  degree  4  with  coefficients  cq) 
and  the  multiplicative  bias  (b  ~  A^(.5,  .5)).  Thus,  we  simulated  values  of  gq  and  b  from  the 
posterior  distribution  (7)  at  each  site,  and  we  obtained  the  following  adjusted  Models-3  values 
(adjusted  value  =  ax  (Models  3) -1-5)  at  the  6  selected  sites:  0.22,  2.90,  1.67,  2.36,  0.96,  and 
1.00.  These  values  are  again  similar  to  CASTNet,  especially  considering  that  the  uncertainty 
about  CASTNet  is  approximately  0.8ppb.  Therefore,  the  methodology  presented  here  for 
combining  spatial  data  is  not  only  useful  for  improving  the  prediction  of  air  quality  but  it 
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is  also  a  tool  for  validating  air  quality  numerical  models  and  for  quantifying  bias  for  bias 
removal. 

Figure  9  shows  the  predicted  values  of  SO2  at  different  locations  in  a  regular  grid,  using 
a  Bayesian  melding  approach  for  prediction  to  combine  CASTNet  and  Models-3  data.  The 
predicted  values  in  Figure  9  are  the  mean  of  the  posterior  predictive  distribution  (eq.  (5)) 
for  the  SO2.  This  graph  looks  similar  to  the  output  of  Models-3  shown  in  Figure  2.  However, 
the  SO2  values  in  Figure  9  are  between  Oppb  and  8ppb,  while  in  Figure  2  the  SO2  values  were 
between  Oppb  and  40ppb.  The  range  of  values  in  Figure  9  is  more  reasonable,  and  it  is  closer 
to  the  range  of  values  for  CASTNet  shown  in  Figure  1.  This  illustrates  the  effectiveness  of 
the  Bayesian  melding  approach  for  correcting  the  bias  in  Models-3.  In  Figure  10  we  have 
6  ensembles  of  SO2  concentrations,  where  each  ensemble  is  a  simulation  from  the  posterior 
predictive  distribution  (eq.  (5))  of  the  SO2.  The  variability  from  ensemble  to  ensemble  is 
due  to  the  uncertainty  in  the  prediction.  Figure  11  shows  the  standard  error  of  the  posterior 
predictive  distribution  at  each  point.  There  is  higher  uncertainty  in  the  Midwest. 

5  Discussion 

In  this  paper  we  have  introduced  a  new  statistical  methodology  for  validation  and  adjustment 
of  numerical  models,  and  for  spatial  prediction  combining  data  with  the  output  from  nu¬ 
merical  models.  In  the  application  presented  in  this  paper,  we  assume  there  is  some  smooth 
underlying  (but  unobserved)  spatial  held  that  measures  the  ’’true”  concentration/hux  of  the 
pollutant  at  each  location.  The  data  collected  at  the  monitoring  sites  are  considered  the 
“true  values”  plus  some  measurement  error.  The  output  of  the  air  quality  numerical  models 
can  be  also  written  in  terms  of  the  true  underlying,  but  unobservable,  process,  with  some 
parameters  that  explain  the  bias  and  microscale  noise  in  the  numerical  models.  The  truth 
is  assumed  to  be  a  smooth  underlying  spatial  process  with  some  parameters  that  explain 
the  large  scale  and  short  scale  dependency  structure  of  the  air  pollutants.  This  spatial  held 
is  represented  locally  as  a  stationary  isotropic  random  held,  but  the  parameters  of  the  sta¬ 
tionary  random  held  are  allowed  to  vary  across  space.  Kernel  functions  are  used  to  ensure 
that  the  held  is  well-dehned  but  also  continuous.  In  this  paper  we  also  take  into  account  the 
change  of  support  problem  that  occurs  when  combining  data  with  different  spatial  resolution. 

Our  objectives  are  model  validation  and  bias  removal  for  the  air  quality  numerical  models, 
and  construction  of  reliable  maps  of  air  pollution  combining  the  output  of  numerical  mod¬ 
els  with  air  pollution  measurements  at  monitoring  sites.  We  validate  air  quality  models  by 
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Figure  8:  The  graph  on  the  left  shows  CASTNet  measurements  for  the  week  starting  July 
11,  1995,  versus  the  modes  and  90%  credible  intervals  of  the  predictive  Bayesian  distribution 
(p(Z|Z,  a  =  0,6  =  1))  given  Models-3  at  the  CASTNet  locations  for  Models-3  validation. 
The  graph  on  the  right  shows  the  CASTNet  measurements  versus  the  modes  and  90%  cred¬ 
ible  intervals  of  the  predictive  Bayesian  distributions  {p{Z\Z^  Z))  derived  from  a  Bayesian 
melding  approach  to  combine  observations  and  Models-3  output.  The  dotted  lines  indicate 
a  90%  credible  region  for  the  CASTNet  values. 
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S02  Concentrations  (Bayesian  melding) 


Figure  9:  Predicted  SO2  concentrations  via  a  Bayesian  melding  approach  to  combine  CAST- 
Net  and  Models-3  data.  This  graph  shows  the  mean  of  the  posterior  predictive  distribution 
for  the  underlying  process  Z  given  CASTNet  and  Models-3. 
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Figure  10:  Simulated  SO2  concentrations  from  the  predictive  distribution  for  the  SO2,  using 
a  Bayesian  melding  approach  for  prediction  to  combine  CASTNet  and  Models-3  data.  The 
simulated  values  are  from  the  posterior  predictive  distribution  for  the  underlying  process  Z 
given  CASTNet  and  Models-3 
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Standard  Error  of  Bayesian  Prediction 


Figure  11:  Standard  error  of  the  posterior  predictive  distribution  for  the  SO2  concentrations 
of  Z,  using  a  Bayesian  melding  approach  for  prediction  to  combine  CASTNet  and  Models-3 
data. 
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obtaining  the  posterior  predictive  distribution  of  the  measurements  at  the  monitoring  sites 
given  the  numerical  models  output.  We  remove  the  bias  the  air  quality  models  by  obtaining 
the  posterior  distribution  of  the  bias  parameters  given  the  measurements  at  the  monitoring 
sites  and  the  numerical  models  output.  We  construct  reliable  maps  of  air  pollutants  simulat¬ 
ing  values  from  the  posterior  predictive  distribution  of  the  true  values  (underlying  process) 
given  the  measurements  at  the  monitoring  sites  and  the  numerical  models  output. 

In  this  paper  we  have  focused  on  air  quality  applications,  but  the  methods  presented  here 
could  be  extended  to  other  problems  where  we  need  to  combine  data  from  various  sources 
and  with  different  spatial/temporal  resolutions.  For  instance,  combining  spatial  data  is  from 
models  and  observations  a  frequent  problem  in  climate  and  weather  prediction. 

Another  approach  to  model  validation  is  to  use  spatio-temporal  models  for  monitoring 
data  to  provide  estimates  of  average  concentrations  over  grid  cells  corresponding  to  model 
prediction  (Dennis  et  al.  (1990),  Sampson  and  Guttorp  (1998)).  This  approach  is  reasonable 
when  the  monitoring  data  are  dense  enough  that  we  can  ht  an  appropriate  spatio-temporal 
model  to  the  data.  In  situations  like  the  one  presented  here,  with  few  and  sparse  data  points 
that  show  a  lack  of  stationary,  the  interpolated  grid  square  averages  would  be  poor  because 
of  the  sparseness  of  the  CASTNet  network,  and  so  treating  them  as  ground  truth  for  model 
validation  would  be  questionable. 

The  spatial  patterns  shown  by  the  air  pollutant  fluxes  and  concentrations  change  with 
location,  so  that  the  underlying  process  Z  with  the  true  values  of  fluxes/concentrations  of  air 
pollution  is  nonstationary  and  standard  methods  of  spatial  modeling  and  interpolation  are 
inadequate.  In  recent  years,  probably  the  most  extensively  studied  method  for  nonstationary 
spatial  processes  is  the  deformation  approach  due  to  Sampson  and  Guttorp  (1992);  see 
also  Guttorp  and  Sampson  (1994),  and  Guttorp,  Meiring  and  Sampson  (1994).  Maximum 
likelihood  versions  of  the  method  were  developed  by  Mardia  and  Goodall  (1993)  and  Smith 
(1996).  In  a  series  of  papers  best  represented  by  Haas  (1995),  T.  Haas  has  proposed  an 
approach  to  nonstationary  spatial  kriging  based  on  moving  windows.  Higdon,  Swall  and 
Kern  (1999)  give  a  model  for  accounting  for  heterogeneity  in  the  spatial  covariance  function 
of  a  spatial  process,  using  a  moving  average  specihcation  of  a  Gaussian  process.  Another 
approach  has  been  developed  by  Nychka  and  Saltzman  (1998)  and  Holland  et  al.  (1999), 
that  extends  the  “empirical  orthogonal  functions”  (EOF)  approach  that  is  popular  among 
atmospheric  scientists.  Here  we  use  a  new  model  for  nonstationary  processes  proposed  by 
Fuentes  (2001a,b),  and  further  developed  by  Fuentes  and  Smith  (2001).  In  this  model  the 
process  is  represented  locally  as  a  stationary  isotropic  random  held,  but  the  parameters  of 
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the  stationary  random  field  are  allowed  to  vary  across  space.  With  this  model  we  are  able 
to  make  inferences  about  the  nonstationary  random  field  with  only  one  realization  of  the 
process. 

The  approach  presented  in  this  paper  gave  us  a  good  understanding  of  the  spatial  struc¬ 
ture  of  the  “true”  concentrations  of  SO2.  This  information  can  be  very  useful  for  designing 
future  data  collection.  Part  of  our  future  work  is  to  use  the  findings  in  this  paper  for  moni¬ 
toring  network  design. 
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